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Abstract. We discuss the effects of open boundary conditions and boundary induced 
drift on condensation phenomena in the pair-factorized steady states transport process, 
a versatile model for stochastic transport with tunable nearest-neighbour interactions. 
Varying the specific type of the boundary implementation as well as the presence of a 
particle drift, we observe phase diagrams that are similar but richer compared to those 
of the simpler zero-range process with open boundary conditions. Tuning our model 
towards zero-range-process-like properties we are able to study boundary induced 
effects in the transition regime from zero-range interactions to short-range interactions. 
We discuss the emerging phase structure where spatially extended condensates can 
be observed at the boundaries as well as in the bulk system and compare it to the 
situation with periodic boundaries, where the dynamics leads to the formation of a 
single condensate in the bulk. 
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1. Introduction 


Stochastic mass transport processes such as the asymmetric simple exclusion process 
(ASEP) or the zero-range process (ZRP) proposed by Spitzer are simple transport 
models for particle hopping aiming to improve the understanding of basic phenomena 
in the dynamics of particles in driven diffusive systems. Generally, these particles 
are abstract and may represent objects from the microscopic to the macroscopic scale 
when combined with appropriate dynamics. It is this relation of abstract particles 
and a multitude of different kinds of dynamics that generates manifold mappings 
to physical processes and phenomena. One such phenomenon that is of particular 
interest to us, is the formation of particle condensates. In fact, dynamics leading 
to steady states in closed, periodic systems where particles form condensates have 
been studied already for the ZRP [T]^ as well as for processes with short-range 
interactions 10 12 . On inhomogeneous structures such as a star graph or scale-free 


networks even the most simple dynamics of uniform hopping can lead to condensation 
at the inhomogeneities 13 -15 . On a homogeneous structure, condensates can emerge 


anywhere in the system as long as the interaction potential falls off sufficiently fast 10 


For a general overview of stochastic transport processes and condensation phenomena 


we refer the reader to the reviews by Schiitz 16 and Evans and Waclaw 20,21 or the 


book by Schadschneider et ah 17 


While the ZRP as well as the extended models can be considered to be driven far 
from equilibrium, their steady state that leads to the condensation remains the same 
as in equilibrium. In fact, in the case of systems with periodic boundaries with particle 
conservation, they are constructed to have this property. This is, however, not a general 
property of transport processes, as can be seen in the exclusion model of Katz, Lebowitz 


and Spohn 18,19 where the stationary distribution may or may not depend on the 
external held depending on the interaction parameters of the model. It is, however, 
also of interest to understand the changes to the condensation process when this steady 
state is broken by replacing the periodic with open boundaries through which particles 
can enter or leave the system, thereby creating a current. In general, this external drive 
and current can lead to phase separation 22 . In fact, for the ZRP, a specihc study has 


been performed by Levine et ah 23 , were among other results phase separation due to 


the introduced boundary drive has been observed. In this paper we seek to extend this 
approach to a stochastic transport process with short-range interactions that feature 
spatially extended condensates in its steady state. This is of interest to us because, in 
contrast to the ZRP, such an extended process is able to interact with the boundary due 
to its non-zero interaction range. As a consequence we are forced to discuss different 
types of open boundaries to grasp their effects on possible condensate formation and 
dynamic phases. Also, instead of using a simpler transport process with short-range 
interactions such as proposed by Evans et ah [^, we decided to employ a tunable 
model 11,12 that can be parameterized to resemble the condensation properties of 
the ZRP as well as extended condensates such as those considered in Ref. fTol. This 
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allows us to compare properties of this model to those of the ZRP discussed in 23 
before going into detail with different types of boundaries. Because the short-range 
interactions in that class of transport processes are strongly related to the fact that the 
steady state of a closed system factorizes over pairs of adjacent sites, we will sometimes 
use the term pair-factorized steady states (PFSS) models although with open boundaries 
a steady state does not necessarily exist. In a previous short note 24 , we have already 
briefly discussed emerging phases and effects caused by the driven open boundaries. 
We did, however, consider only one specific type of open boundaries and were severely 


limited by the employed numerical method. In a recent short communication 25 , we 
sketched an improved simulation setup and discussed for this special case the phase 
diagram and transition dynamics between the phases in more detail. In particular, we 
pointed out that not only the details but, in fact, the very existence of phases depends 
on the choice of interaction with the boundary. We therefore would like to complete 
the picture with that versatile numerical approach and an emphasis on the point that 
the specific interaction details at the boundaries have significant impact on the system’s 
phase diagram. 

The remainder of this paper is organized as follows. In the next section we will 
briefly introduce the zero-range process as well as the tunable short-range interaction 
stochastic transport model and define the considered types of open boundaries. In the 
third section we describe the used numerical methods and motivate our choice for a 
kinetic Monte Carlo algorithm. In the fourth section we will discuss our results, first 
making a comparison with the zero-range process and then discussing emerging phases 
and properties in detail with short-range interactions turned on. Finally, we summarize 
our findings in the fifth section. 


2. Stochastic transport processes with open boundaries 

The basic particle-hopping stochastic transport process consists of a one-dimensional 
lattice with L sites and a gas of M indistinguishable particles. Each site i of the lattice 
can be occupied by any number rrii = 0,..., M of particles, where M = 
the total number of particles. These particles can leave their sites i with a rate Ui and 
then jump to one of the adjacent sites. A target site is randomly selected among the 
neighbours with respect to the strength of asymmetric hopping given by the parameters 
p and q so that the actual rates of particles hopping to the sites right (i -|- 1) and left 
(i — 1) of the departure site become pui and qui, respectively, as indicated in Fig. At 
the boundaries, particles enter or leave the system. We define exchange rate parameters 
Pin, Qin and Pout, Qout to control particle currents into and out of the system. The exact 
mechanisms of injection and removal of particles at the boundaries are discussed further 
below where we define the specific properties and implementations of the boundaries. 
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Figure 1. Schematic representation of the dynamics of a particle hopping process on 
a one-dimensional lattice with L sites, hopping rate Ui and drift parameters p, q. At 
the boundary sites, the drift parameters towards the boundary are replaced with the 
removal parameters Pout and gout as indicated. Likewise, the rate of particle injection 
at the boundaries is given by the parameters pin and gin, respectively. 


2.1. Zero-range proeess 

The hopping rate function Ui determines the dynamics of the particles in the system. 
For zero-range processes it must only depend on the occupation number of the departure 
site which results in a local-only interaction term. An example for a specific ZRP, that 
is of particular interest in the context of this work, is the condensation model with 
hopping rates 

u{m) = 1 b/m, (1) 

where a single-site particle condensate spontaneously emerges for 6 > 2 in the steady 
state of the periodic system when the particle density exceeds the critical density 
Pcrit,ZRP = ^/{b — 2) }2^ . 


2.2. PFSS process 


In this paper, we consider a more generic model where the hopping rate function depends 


also on the number of particles on the adjacent sites as proposed by Evans et ah 10 


With an appropriate choice of the hopping rate function, spatially extended condensates 
emerge in the periodic system due to the nearest-neighbor interaction. One generic 
choice of the hopping rate function of this process reads 

g{mi - l,mj) 


Ui = Ylu{mi,mj) = JJ 
(bi> {i,j} 


g{mi,mj 


( 2 ) 


making the interaction potentials between adjacent sites {i,j) sites isotropic for 
symmetric weight functions g{m,n) = g{n,m). By construction this results in a pair- 
factorized steady state probability distribution of the form 




^M,L n ^ rtii,M 1 

irj) 


( 3 ) 


as long as the number of particles is conserved. Here {m} gives a complete state, Zm,l 
normalizes the steady state similar as the partition function in an equilibrium system 
and the Kronecker symbol fixes the particle number. 
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To be able to compare to the work on ZRP with open boundaries 23 as well as 


to consider a process with effective long-range interactions, we use the hopping rate 
function generated by the tunable interaction terms 


g{m, n) = exp 


m 


n 


/3 


-im? + m?) 


(4) 


proposed by Waclaw et ah 11,0^. The weights g{m,n) consist of a zero-range 


interaction term tuned by the parameter 7 and a nearest-neighbour interaction term 
tuned by the parameter (3. The steady state of this process features the formation of 
particle condensates of various properties depending on the values of these parameters. 
A critical density and thus condensation phenomena exists for 0 < 7 < 1. The 
condensate then assumes one out of three qualitatively distinct forms that strongly 
influences the model’s dynamics: a single-site peak for /? < 7 , an extended rectangular 
shape for 7 < /3 < 1 or a smooth parabolic shape for (3 > 1 11 12 27 . Most important to 


our purpose is the ability to reproduce single-site condensates similar to those observed 
for the ZRP with hopping rates Q for /3 < 7 and 7 < 1, as well as spatially extended 
smooth condensates for (3 > 1 and 7 < 1 similar to those observed in Ref. 10 . 


2.3. Open boundaries: Mechanisms of particle exchange and external drive 

For the zero-range process the implementation of open boundaries is straightforward 
because there is no interaction other than particle exchange. For the considered model 
(§-(§, on the other hand, the type of interaction at the boundary sites i = 1 and i = L 
has to be chosen explicitly due to its non-zero interaction range. We will consider and 
discuss two main types of implementations. 

Our hrst approach is to interpret the system as isolated and discard the interaction 
terms of the bonds that cross the boundary as follows from the factorization of weights 
for arbitrary graphs in Eq. ([^. In the following, we will refer to this type as loose 
boundaries. As a second approach we consider the system to be embedded in a larger 
system with a separation that hinders particle movement like a membrane. Here, we do 
not discard the interaction term for the bond that crosses the boundary as it reflects 
the interaction with some mean-field occupation outside of the considered system by 
setting the external particle occupation to a constant value m^o. In contrast to the 
first approach, we refer to this type with the term fixed boundaries. This results in the 
hopping rates 


n(mi,m 2 ), u{mL,rnL-i) for/oose boundaries, 

Ml, Ml = (5) 

u{mi,m 2 )u{mi,moo),u{mL,mL_i)u{mL,moo) ior fixed boundaries. 


Additionally we consider two types of particle removal at the boundary sites that 
differ in the way the hopping rate determines the rate of particles leaving through the 
boundary. Intuitively, the rates of removal are Wigout at the first and u^Pout at the 


last site. This mechanism is used for the ZRP in Ref. 23 as well as in our own prior 
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study j^. Because normal hopping is involved for particles leaving the system, we use 
the term hopping removal to refer to this. Here, we also consider a second removal 
mechanism for particles that is symmetric to the mechanism of particles entering the 
system that occurs at a constant rate. For the latter mechanism we use the term 
constant removal. Because at the boundary sites the rate of particles leaving the system 
is decoupled from the hopping rate ui and UL,'we replace them with ul and define 

* * I Ml, Ml forremoval, 

Ui,Ul = < ( 6 ) 

1 ^1, 1 for constant removal, 

so that the removal rates become M^^out for the first and M^Pout for the last site. Rates of 
particles moving towards the bulk system remain unchanged {uip and ulQ, respectively). 

In the following, we will only consider exchange rates that in general reflect the 
external drive of the system. That is, for symmetric dynamics {p = q = 1/2) the 
exchange parameters are identical at both boundaries (pm = q'in,Pout = Qout), while for 
totally asymmetric hopping (p = l,g = 0) the exchange is restricted to that spatial 
direction as well (gin = gout = 0). 


3. Numerical simulation methods 


The usual approach to simulating the dynamics of a stochastic transport process as a 
Markov chain is as follows: first propose a random departure site i, second compute 
the acceptance probability for the hop from the hopping rate u* and third decide 
whether a particle hops to a randomly choosen neighbour. To compute an acceptance 
probability it is required, however, that the hopping rate function can be normalized for 
any permissible local combination of occupation numbers. This normalization basically 
results in a change of the simulation time scale by the normalization factor. 

While this is possible both for hopping rate functions with an upper bound, such 
as Eq. ([^ and those proposed by Evans et ah (^, or when a maximum rate is known 
due to conservation of the total number of particles M{t) = const, it becomes inefficient 
for increasingly separated intrinsic time scales of slow and fast events and thus results 
in a large ratio of rejected updates. 

The hopping rates of the model considered here, however, do not have an upper 
bound in the regime f3 > 1. In fact, the required normalization constant would grow 
roughly as the square of the number of particles in the system, so that the approach 
to directly simulate the dynamics as sketched above cannot be used here. To work 
around this limitation as well as to improve efficiency in the presence of fast and slow 
events we employ a rejection free kinetic Monte Carlo (KMC) algorithm introduced 
as the direct method by Gillespie 28,29 for the simulation of coupled rate equations. 


While the method was designed for small chemical systems with few reactions, it can 
be made fit for efficient simulation of larger systems with some optimizations. The idea 
of the method is similar to those of other rejection-free KMC methods such as the n- 


fold way algorithm 30,31 in that an update consists of selecting an event according to 
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its specific rate Tk relative to the total rate V = of all possible events, execute 

it and update the system time t —)■ t + At, where At is the waiting time to generate 
this event. Finally the list of possible events and their assigned rates F*, are updated to 
reflect the new state of the system. In the direct method, the event k is picked by the 
relation < Fxi < and the exponentially distributed time increment 

At is determined as At = — (Inx 2 )/F using two uniformly distributed random numbers 
Xi^X 2 € [0,1). For our model these events are all the particle transfers between adjacent 
sites (i, j) with their rates Uip, Uiq for 1 < i < L, uip, ulQ and the removal and injection 
of particles at the boundaries with rates Wi^out, WiPout and Pm, fe, respectively. 

The search for the appropriate event is easily improved by using a binary search 
combined with multiple levels of search, where events that originate from 


tree 32 


the same site are grouped in the first search level, and then resolving to one specific 
event for that site. The step to update the rates of events is efficiently implemented 
by taking into account which events’ rates actually need updating based on which 
neighbouring sites were involved in the last step of the Markov chain. This is basically 


an application of the proposed update principle of the next-reaction method 33 which 


involves building dependency graphs between events and rate recalculation to achieve 
this. As an additional advantage to simpler methods, the time scale of simulations 
becomes equivalent to the physical time scale, thus making it obsolete to define artificial 
time scales in terms of sweeps or local updates. 

Because of the continuous time simulation method, we cannot compare CPU time 
per full update, but per unit physical time of the simulated system. For a lattice size of 
L = 256 sites, the simple simulation for the ZRP takes Sips, somewhat faster than the 
KMC method with 37ps for one unit of model time. In a situation with slow and fast 
events, however, the simple method becomes slower proportional to the ratio of large to 
small rates, while the performance of the KMC method does not suffer from that. In 
our simulations the typical inprovement factor with only M 100 would be around 50 
but growing roughly with the square of the total particle number. 

To compute the observables for the phase diagrams, we simulated at least 25 replicas 
for each point (pin,Pout), and between 50 and 100 for points near the transition lines. 

4. Results 


4-1. Open boundary effects in the zero-range process like regime 

We start the discussion of boundary drive induced dynamical phases with a look at 
the zero-range process with hopping rates Q. Most notably, for 6 > 2, it features 
spontaneous symmetry breaking and the formation of a single-site particle condensate 
in its steady state for periodic boundaries. That is, a single site contains a finite fraction 
1 — Pcrit/p of all particles, where pcnt = — 2) is the critical density that is assumed 

on average in the rest of the system. Effects of open, driven boundaries on this model 


have been studied and discussed by Levine et ah 23 . For the ZRP, we will use the 
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Figure 2. Scaling parameter a of the total number of particles in the ZRP-like 
regime for /3 = 0.4,7 = 0.6 and the ZRP for 6 = 5. Low values of a « 0 indicate 
a stable, fluctuating total number of particles while a = 1 shows linear growth in 
time. Diagrams (a), (b) correspond to the ZRP-like model, (c),(d) to the ZRP with 
symmetric and totally asymmetric dynamics respectively. 


parameter 6 = 5 that results in a critical density of pcnt = 1/3. In the following we 
tune the coupling parameters [3 and 7 of our model Eq. (|^ to a regime where the 
condensation process with periodic boundary conditions has similar properties in the 
steady state as the ZRP and compare some of the properties to those found for the ZRP. 
With the choice of parameters /? = 0.4 and 7 = 0.6 in Eq. Q a single-site condensate 
and critical density of pcnt = 0.302 ± 0.006 similar to that of the ZRP is expected for 
periodic boundary conditions (27|. We use loose boundaries with hopping removal, the 


latter of which is the same as that used for the ZRP in 23 


Two major phases are expected for the ZRP with open boundaries 23 . First, a 


steady state with a thin homogenous particle gas, and second, a phase with aggregate 
condensates formed at one or both boundaries that act as particle reservoirs and can 
influence the bulk system in between. To distinguish between these phases we measure 
time series of the total number of particles M{t) and the bulk density Pbuik- It is useful 
to determine a scaling exponent a for M{t) assuming it roughly follows a power-law 
M{t) ~ to hnd the two phases. The bulk density is estimated as the average particle 
density in the bulk of the system 


*bl 

Pbuik = w—V-— rUi, where m* > 0 V i < ibf and i > ibi, (7) 

1 + *bi - ^bf ^ 

*=*bf 

that is, the region from the hrst (ibf) to the last (ibi) unoccupied site in the system. 

From the plots of the scaling exponent a given in Fig. |^as well as the bulk density 
Pbuik given in Fig. [^we can clearly identify the same phases for this regime of our tunable 
model Q as of the ZRP. In both models there is a particle gas phase (G) with a low 
stationary particle density p = pbuik that increases towards the transition line to the 
aggregate condensate phase (A). There, large numbers of particles aggregate at the hrst 
and/or last sites of the system. While the latter phase is homogeneous in systems with 
symmetric hopping (p = q = 1 / 2 ), sub-phases Ain, Aout and A, where the aggregate 

































9 


PFSS = 0.4,7 = 0.6) 




Pbulki- I I I Pent = 0.302 ± 0.006 

0 0.1 0.2 0.3 0.4 


ZRP (6 = 5) 




/^bulk P I I Perit — 1/^ 

0 0.1 0.2 0.3 0.4 


Figure 3. Bulk density pbuik in the ZRP-like regime for j5 = 0.4 ,7 = 0.6 and the ZRP 
for 6 = 5. The critical densities for the PFSS and ZRP are Pcrit.PFSS = 0.302 ± 0.006 
and Pcrit.ZRP = 1/3, respectively. Except for the non-criticality of the bulk system, the 
dynamical phase diagram of the tunable system is to a high degree similar to that of 
the ZRP. 


condensate forms at the first site t = 1, the last site i = L, or both, can be identihed 
for asymmetric hopping (p 7 ^ g), see Figs. (b), (d), and (b), (d). The aggregate 
condensates at the boundary sites act as reservoirs for particles entering (Ain) and 
leaving (Aout) the system, effectively regulating particle flux through these sites. This 
can be seen well for totally asymmetric hopping, where in A^ the bulk density assumes 
the value Pbuik = 0.15 ± 0.04 < pent in the tunable model and Pbuik ~ Pent = 1/3 
for the ZRP. In Aout the reservoir cannot act on the bulk system, so that the bulk 
system is still a particle gas. For symmetric hopping, however, both phases combine 
and aggregate condensates at both boundaries act on the bulk system, increasing its 
density to pbuik = 0.20 ± 0.04 for the tunable model, still below criticality. The bulk 
system of the ZRP remains critical and long-lived bulk condensates sometimes emerge, 
which results in large values of the bulk density as shown in Fig. |^c). 

To understand the formation of a condensate at the influx boundary for totally 
asymmetric hopping a simple biased random walk in the occupation number of the hrst 
site can be considered 23 . Because the drift pin — (1 -|- b/m) of the walker becomes 
positive for influx rates Pin > 1 and sufficiently high occupation of the hrst site, a stable 
condensate can emerge at that site. In fact, we observe a power-law dependence of the 
waiting time on the inhux rate Pin until the Aju-condensate forms after a quench to the 
aggregate condensate phase. This corresponds to the hrst-passage time of that random 
walk process to a sufficiently high occupation number where it has positive drift. For 
symmetric hopping, the argument is similar but results in a diagonal transition line to 
the aggregate condensate phases for pin > Pout because particles may leave directly after 
entering. For a more detailed discussion of this argument, also with respect to partially 
asymmetric hopping for the ZRP, we refer to the original work of Levine et ah |^ . 

The same argument can be applied to our model ([^-(|^ in the regime // < 1, when 
the weak short-range interactions are negligible. The resulting hopping rate at the hrst 
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site becomes 


Ml = exp 



(mi — l)"^) + |mi|^ — |mi — 1|^ 


approximated by 


Ml = exp 


^'yrnj ^ + (3m^ ^ 


( 8 ) 

(9) 


for large values of mi. This in turn approaches the value 1 for large occupation numbers 
mi as long as /3 < 1, that is for the entire single-site and rectangular condensate regimes 
of the model. The drift of the first site’s occupation becomes positive for the same value 
of Pin > 1 and yields therefore the same transition line as for the ZRP. 

The formation of the aggregate condensate at the outflux boundary at site L is 
easily understood in the totally asymmetric case with a similar argument as above. For 
Pin < 1 all particles eventually reach the last site L. If the removal rate is smaller than 
the rate of particles arriving at the site pout < Pin, the drift of the occupation number 
rriL becomes positive and an aggregate condensate emerges. 


4-2. Open boundary effects in the extended condensate regime 


The goal of this section is to identify the qualitative phase structure of the described 
transport model depending on the strengths of particle exchange at the boundaries 
with respect to the considered types of boundary conditions. Within this section, the 
interaction parameters of the tunable transport process Eq. (|^ are fixed at /? = 1.2 and 
7 = 0.6, setting it into the regime of smooth parabolic condensate shapes of the periodic 
system 27 . The critical density for these parameters in a comparable system with 
(L = 256, M L) is pcrit ~ 0.3 due to finite-size effects and decreases to 0.125 ± 0.009 
in the limit of large systems. 

Based on the phase structure of the ZRP given in Ref. 


23 and numerically 


reproduced for the ZRP and a short-range interaction transport model with smooth 


condensates in this and our own previous work 24 we expect to some extent a similar 


phase diagram. Therefore, to identify the phases, we continue to use the time series of 
the total number of particles M{t), its scaling exponent a and the bulk system particle 
density pbuik introduced in the previous section. An example of the total mass versus 
time M{t) ~ for loose boundaries and constant removal along with numerically 
determined values of a is given in Fig. Additionally to these quantities we record the 
microstates of the systems at regular intervals, so that we can compute other quantities 
such as the occupation number profiles that we use later. 

As shown in Fig. the scaling exponent a identifies regions with distinct values of 
a 0 and a ~ 1, that is, stationary as well as linear growing total numbers of particles 
M{t), respectively. For constant particle removal, additionally the value a ~ 0.6 is 
observed on the transition line between these former regions. Together with the data 
for the bulk density shown in Fig. we are able to identify candidates for gas phases 
with low values of a = 0 and Pbuik, and aggregate condensate phases where a = 1 and 
low values of pbuik are observed. Additionally a phase with stationary particle count 
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Figure 4. Average total number of particles M{t) ~ t“ for loose boundaries 
with constant removal of particles and symmetric dynamics (pin = 9in;Pout = 9out) 
determined from 25 replicas. There are four distinct groups of curves: linear growth 
of M (t) for Pin > Pout, approximate square root growth for pin = Pout and two groups 
with stationary particle numbers, both for pin < Pout- The straight grey lines indicate 
the different observed types of scaling. The scaling parameter a as determined from 
the average slope in the log-log-plot in the interval 10^ < t < 10® is given right to the 
respective key symbol. 


but relatively large bulk density is found in between those for constant removal and 
symmetric hopping. To exactly identify the type of phase a system is in at any given 
parameterization (pin,Pout) we use graphical representations of the individual systems’ 
evolution of microstates over time such as shown in Fig. and averaged occupation 
number prohles computed from many individual trajectories shown in Fig. Combining 
this information we are fully set up to identify the regions in the phase diagrams and 
discuss their properties in the following subsections. 

4 . 2 . 1 . Particle gas phase (G): For the considered types of interactions at the 
boundaries, a particle gas phase (G) as observed for the ZRP and ZRP-like regime of 
the tunable model exists. Likewise it features a thin gas of particles hlling the complete 
system. It is observed for small enough values of influx rates p\n und large enough 
outflux rates for symmetric hopping pout, so that particles can directly enter and leave 
the bulk system. In the gas phase the system can be thought of as being part of a larger 
periodic system. The stationary particle density p = pbuik (where a = 0 as shown in 
Fig. 1^ increases with stronger drive at the boundaries towards the critical density of 
the steady state system as shown in Fig. For the small system sizes we considered 
to obtain most of our data, the critical density is signihcantly influenced by hnite-size 
effects. That is, for a small total number of particles as observed in the gas phase, the 
critical density pcut ~ 0.3, were condensation hrst occurs, is considerably larger than the 
large system limit pcnt = 0.125±0.009 which is approached with increasing total particle 
number as shown in Fig. |^a) for the ZRP and ZRP-like model and Fig. [^b) for the 
extended condensate regime of the short-range model with /? = 1.2,7 = 0.6. Note that 
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Figure 5. Average scaling exponent a of the total number of particles in the 
system and phase boundaries for the various types of boundary conditions. The 
boundary type is given by combination of the labels on the left and top margins, 
e.g., (f) loose boundaries with constant removal and totally asymmetric dynamics 
{p = l,q = qin = 9out = 0). The additional spanning condensate (SC) phase, which 
features a single stationary bulk condensate of maximal width, as well as the spanning 
fluid (SF) phase, where the system absorbs new particles, in panels (e) and (g) will be 
discussed further below in the text. 


the observed bulk density pbuik for the smaller systems is above the asymptotic value of 
Pcrit, but still below pcut of fhe hnite system as it should be. The values for the critical 
density given in Fig. [^a) were determined as the background density of a periodic 
system with overall particle density signihcantly above the condensation threshold. For 
a very similar ZRP with hopping rates u{m) = 1 + b/mP, such hnite-size effects have 


already been observed 34,35 


4.2.2. Aggregate condensate phases (A): For sufficiently large influx rates Pin, particles 
tend to become adsorbed at the boundary and form aggregate condensates. This is 
observed for any of the boundary types. As in the ZRP and the ZRP-like regime, 
this phase consists actually of three regions with aggregate condensates at the influx 
boundary, the outflux boundary and at both boundaries that exist individually for 
totally asymmetric hopping and mix to a single uniform region for symmetric hopping. 

An example time series of an inbound aggregate condensate absorbing entering 
particles is shown in Fig. [^b). The bulk density Pbuik as observed in Fig. is 
consistently below the asymptotic value of the critical density as determined in Fig. [^b). 
The aggregate condensates show individual shapes depending on whether they absorb 
inbound or outbound particles as well as on the type of boundary {loose, fixed) and 
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Figure 6. Particle density pbuik in the bulk system for given boundary drives and the 
four considered types of open boundaries. The plot values are cut off at pbuik = 0.25 
to retain readability for the system with constant particle removal, where due to the 
large bulk condensate the density is increased by orders of magnitude. In the dotted 
regions, the bulk density remains undetermined as the bulk condensate has been in 
contact with the system boundaries in all simulated replicas. The hatched region in 
panels (b), (f) and (h) marks a transition region between the Aout and A phases. The 
individual plots for each boundary type are labeled as in Fig. 
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Figure 7. Example time series: (a) Gas phase (G, loose/hopping, p = q = 1/2, 
Pin = 1.25,Pout = 1), (b) boundary aggregate condensate (Ain, fixed/hopping, p = 
l,q = 0,pin = Pout = 2), (c) spanning bulk condensate (SC, loose/constant, p = 
9,Pin = 1-25,Pout = 1-5), (d) intermediate spanning fluid phase (SF, loose/constant^ 
P = 9,Pin = Pout = 2). 
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Figure 8. Average occupation number profiles rrii at simulation time t = 10® for the 
considered types of boundary conditions. Plot symbols refer to influx rates, colors to 
outflux. The boundary types are from top (a, b) loose with hopping removal, (c, d) 
fixed with hopping removal, (e, f) loose with constant removal and (g, h) fixed with 
constant removal. The left-hand side plots represent results for symmetric dynamics 
{p = q = 1/2,Pin = (7in,Pout = 9out), the right-hand side results for totally asymmetric 
dynamics {p = l,q = q-m = qout = 0). To improve readability, not all points are plotted 
as symbols. To compute these profiles, 25 to 40 trajectories were used. 
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Figure 9. Finite-size effects of the critical density for low to high overall density 
in a closed periodic system of different sizes of L = 256,512,1024 and 2048 sites for 

(a) the ZRP and the effectively ZRP-like process with with /3 = 0.4,7 = 0.6 and 

(b) short-range interactions (3 = 1 . 2,7 = 0-6- For sufficiently large particle numbers 
M beyond the visible local maxima, condensates emerge and the bulk density pbuik 
becomes the critical density Pcrit- Points represent data obtained from simulation of 
the steady state (10® Monte Carlo sweeps), lines show the fitted finite-size scaling law 
Pbuik = Pcrit+oM“^ with b « 0.40, where the actual critical density pcrit = 0.125±0.009 
is approached for increasing system size L and particle number M. 


particle removal {hopping, constant). The qualitative shape of inbound aggregate 
condensates Am [Fig. ib. d, f, h)] closely resembles the steady state condensate 
shape of the model with periodic boundaries. The shape of the outbound condensate 
Aout [Fig. gf, h)] has a relatively steep increase of occupation numbers towards 
the boundaries but becomes almost flat when approaching the transition to A with 
increasing influx pin. In this transition zone [hatched area in phase diagrams of 
Fig.gb, f, h)] the aggregate condensates show significantly increased widths with respect 
to mass, so that merging of both aggregate condensates is observed very early compared 
to the region A. As a low density bulk cannot exist, these transition regions are also 
dotted in the phase diagrams of Fig. 

The difference in condensate shapes caused by the boundary types can be seen by 
comparing the respective profiles row by row in Fig. |^and in fact is visible also in the 
large bulk condensate phase where it approaches the boundary as discussed below. For 
loose boundaries, the prohle starts and ends at high occupation numbers with zero slope 
at the boundary. Since there is no interaction beyond the boundary, its shape is very 
close to one half of a steady state shape. With fixed boundaries, the condensate shape is 
forced to lower occupation numbers towards the boundaries by the interaction term with 
the mean-field occupation moo = 0. Also, the maximum occupation of the condensates 
becomes lower which leads to increased condensate widths because the total mass of the 
condensates is comparably independent of loose or fixed boundaries. The mechanism 
of particle removal at the boundaries seems to affect the aggregate condensates only to 
the extent that their rates of growth are changed. Constant removal lets more particles 
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(a) 



(b) 


Figure 10. (a) Ratio of replicas with an emerged aggregate condensate after a waiting 
time tq^a after a quench to the aggregate condensate phase (101 for influx rate values 
of 1.4 < Pin < 1.6 and constant outflux rate Pout = 1-5 for totally asymmetric hopping. 
For this plot, from N = 200 up to 800 replicas per influx value were used with system 
size L = 256. Larger systems give qualitatively the same result, (b) Log-log plot of 
the excess influx rate pin — Pin.crit versus the half-value waiting time ti/ 2 , where half 
of the replicas have developed an aggregate condensate. Symbols represent numerical 
data determined from simulations of = 200 system replicas of size L = 256, lines 
show the fitted scaling law (111. 


escape the system and therefore results in much smaller aggregate condensates. 

To estimate the slope and position of the transition line to the phase A, we perform 
quenches of a large number of replicas of systems to several values of the influx rate 
Pin > Pin,crit beyoud the transition line, where Pin,crit is the critical influx rate for the 
given parameterization. For any value of the “depth” pin — Pin,crit of the quench into the 
phase we then measure at several times the transition ratio of the gas phase 

1 ^ 

-Ptrans(l'G—> a) E H {M{t) - Mthresh) , (10) 

i=l 


where H{x) is the Heaviside function, N is the number of replicas for a given value of 
Pin and Mthresh = oPcrit-h is a threshold mass (with a > 1) to detect the transition to the 
aggregate condensate phase. This transition ratio is related to the survival probability 
of the gas phase as Psurv('r) = 1 — Ptnmsij)- The determined values for loose boundaries 
and hopping removal are given in Fig. [To|(a). The power-law scaling of the transition 
time becomes evident when looking at a hxed transition ratio Ptrans ~ 0.5 in Fig. [Io|(a): 
approximately equidistant increases of the depth of quench halve the mean waiting time 
to reach that transition ratio. This is illustrated in Fig. 10 b), which shows the scaling 
of the depth of quench versus the waiting time for transition of 1/2 of the replicas. The 
values determined for this scaling relation 


Pin - Pin.crit /I (H) 

near the transition line as well as the critical influx rates Pin,crit sxe given in Table 
Additionally to the considered tunable model we also considered the ZRP to check our 
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Table 1. Numerically determined values for the scaling relation ( [Tl] ) of the quench 
depth Pin — Pin.crit versus the waiting time ti /2 to transition for 1/2 of replicas. These 
values are determined for totally asymmetric hopping and hxed outflux rate Pout = 1-5, 
where the critical influx rate does not depend on Pout- 




scaling exponent k 

critical influx rate Pm.crit 


1.2 ,7 = 0.6, fixed 

0.210 ±0.064 

1.426 ±0.018 


1.2 ,7 = 0.6, loose 

0.200 ±0.027 

1.459 ±0.013 


0.4 ,7 = 0.6, loose 

0.233 ±0.031 

1.0337 ±0.0038 

ZRP, 6 = 5 

0.218 ±0.014 

1.0055 ± 0.0024 


methods. As for the PFSS, we considered the waiting time nntil a condensate attached to 
the inflnx bonndary neglecting formation of droplets in the bnlk system. We only rarely 
observed the case were a droplet forms away from the bonndary and grew to become 
the aggregate condensate. We did not observe the formation of stable condensates in 
the bnlk. We wonld like to snggest that the sitnation for bnlk condensate formation 
here is indeed different to that of the ZRP with periodic bonndaries where coarsening 
sets in immediately (c.f. Ref. [^). The phase transition conld appear, when a droplet 
in the bnlk grows fast enongh to become immobile and attach to the bonndary instead 
of diffnsing or leaving the system. From onr observations, the aggregate condensates in 
fact formed at the bonndary. 

Using these observations of the scaling of transition times to the aggregate 
condensate phase, we reprodnced the valne of the critical inflnx rate Pm.crit = 1 for 
the ZRP (see Ref. j^) as well as determined the scaling exponent k, c.f. Table 

We would also like to point out that the exponents observed for the scaling relation 
of the “depth” of quench with the transition time are within their statistical errors 
identical {k ~ 0.22 ± 0.02), although different models (tunable PFSS and ZRP), 
couplings {fixed and loose) or interaction at the boundary {hopping and constant 
removal) are considered. The physical meaning of this scaling exponent for a quench 
from the gas to the aggregate condensate phase is the connection of the “depth” of 
the quench into the new phase to the time it takes until it manifests in the system (or 
the survival time of the old phase). The value of the scaling exponents here hints at a 
universality of this transition. Possibly related to this is the global persistence scaling 
exponent 6, which describes the distribution of survival times P{t) ~ of the old 


phase after a quench to a critical point 36 38 in non-equilibrium systems. However, 
we would like to postpone this interesting question to future work as we could not yet 
address it properly. 

The nature of the aggregate condensate phase is different to the regime with 
negligible short-range interactions (e.g., fi = 0.4,7 = 0.6), however. This becomes 
clear when looking at the argument of the hnite biased random walk of the occupation 
of the boundary sites. The hopping rate at the hrst site with an unoccupied neighbor 
does not decrease or even approach a stationary value when the occupation number mi 
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increases so that a single site aggregate condensate cannot form. With a neighbor of 
similar occupation, however, the hopping rate Ui has a local minimum for a non-zero 
occupation number, so that for this configuration of the first sites there is a positive drift 
of the occupation number in the random walk argument. Due to this interaction with 
the sites in the system, a spatially extended condensate aggregates at the boundary. 

4.2.3. Spanning condensate phase (SC): With symmetric hopping and constant 
particle removal, an additional phase featuring a single large condensate emerges 
intermediate between the gas and aggregate condensate phases. The condensate spans 
the bulk of the system almost approaching the boundaries. As in the gas phase, the 
total number of particles M{t) and therefore the condensate mass is stationary. This is 
already visible from the high value for the stationary total number of particles given in 
Fig. I^for rates {pm,Pout) = (1-0,1.5) and (1.5,2.0). An example time series leading to 
such a spanning bulk condensate is shown in Fig. [^c). The resulting average occupation 
profile for large times is shown in Fig. [^e, g). There it is visible that towards the phase 
boundary to the aggregate condensate phase the bulk condensate starts to touch the 
boundary sites. This becomes clear in Fig. [^e, g) where the dotted region indicates that 
the measurement of the bulk density as defined in Eq. ([^ cannot be achieved because 
the bulk does not exist. 

To understand why this additional phase occurs with the constant removal 
mechanism, we suggest that for single particles at the boundary as it ocurs in the 
gas phase, the removal rate is much lower than with the hopping removal mechanism. 
With asymmetric hopping this leads to emergence of an outward boundary aggregate 
condensate (Aout) as evident from Fig. it. h). For symmetric hopping, however, this 
leads to an increase of the bulk density above the critical density and thus the formation 
of a bulk condensate. This condensate then absorbs particles until its stable maximum 
size is reached. When the particle influx rate pin becomes large enough to create 
aggregate condensates, the bulk condensate connects to the boundaries resulting in a 
flat occupation profile. We mark this transition with a zig-zag line denoted as spanning 
fluid (SF) phase in Figs. and The total number of particles in this transitionary 
state grows roughly with a ~ 0.6 as shown in Figs. |^and|^ 


5. Conclusion 


We systematically studied the effects of open boundaries and external drive in a 


stochastic transport process with tunable short-range interactions 11, 12 far from 


equilibrium. To do so in a meaningful and systematic way we proposed four different 
boundary types distinguished by the type of interaction and the mechanism of particle 
removal at the boundary. The interaction at loose or fixed boundaries reflects the 
non-existence or existence of an interaction term with a mean-field occupation across 
the boundary site. For leaving particles, additional to hopping removal, we propose 
constant removal for reasons of the symmetry of particle exchange. We considered 
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the four types of open boundaries generated by combinations of these properties and 
determined the respective phase diagrams for both symmetric and totally asymmetric 
hopping dynamics. 

We successfully applied the direct KMC technique, a continous-time Monte Carlo 
method to the system to control the unbounded values of the hopping rate of the chosen 
dynamics for /I > 1 where in the steady state of the closed (periodic) system extended 
condensates of smooth parabolic shape form. This, however, also mitigates the need to 
artihcially formulate update sweeps that coordinate regular particle hops with particle 
injections and removals because all events can be treated equally only according to their 
rates leading to a uniform time scale. 

For negligible strength of the short-range interaction < 7 ) we found a phase 
diagram that is essentially equivalent to that of the ZRP condensation model ([^ as 
discussed in Ref. 23 . A homogeneous particle gas phase and an aggregate condensate 


phase make up the phase diagram. The transition mechanism between these can be 
understood the same way as for the ZRP. 

When the short-range interactions become important (/3 > 1) we hnd an enriched 
phase structure. The particle gas phase is identical to that in the prior models. In the 
aggregate condensate phases, however, spatially extended condensates emerge at the 
boundary sites with envelope shapes that adapt to the predominant flux of particles 
in or out of the system in case of asymmetric dynamics. The interaction at loose or 
fixed boundaries, while not changing the phase structure qualitatively, does have a 
signihcant effect on the transition lines between phases as well as the properties in the 
aggregate condensate phases. In the case of fixed boundaries, this is very obvious in the 
deformation of the aggregate condensates at the boundary sites. With the constant rate 
particle removal mechanism, however, we observed the emergence of a new intermediate 
phase featuring a dominant bulk condensate between the particle gas and aggregate 
condensate phases. To obtain a precise value for the critical influx rate that separates 
phases G and Ain for totally asymmetric hopping, we analyzed survival times of the gas 
state for different quenches to A^. An interesting observation in this analysis was the 
identity of the scaling exponents of the relation between distance to the transition line 
and half-value survival time across different models, coupling strengths and considered 
boundary types. While we could not yet identify the cause of this scaling, it appears to 
be a universal property for this type of phase transition. As a future project it would 
be interesting to work out its possible relation to the global persistence scaling. 
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